Chapter 4 Diversity analysis

load("data/data.Rdata")
load("data/gifts.Rdata")
load("data/beta_27022025.Rdata")

genome_counts_filt <- genome_counts_filt %>% 
  column_to_rownames(var = "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>% 
  select(where(~ sum(.) > 0)) %>% 
  select(-AC85) %>% 
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)

table <- genome_counts_filt %>%
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  rownames_to_column(., "sample")
  
master_table <- sample_metadata %>%
  mutate(sample=Tube_code) %>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>% 
  filter(Tube_code %in% table$sample) %>%
  mutate(treatment = sub("^\\d+_", "", time_point))%>% 
  left_join(., table, by=join_by("Tube_code" =="sample"))

sample_metadata <- master_table %>% 
  select(1:13)

genome_counts_filt <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))%>% 
  rownames_to_column(., "genome")

genome_counts_filtering <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
#match_data(genome_counts_filtering,genome_tree)

genome_metadata <- genome_metadata %>% 
  filter(genome %in% genome_counts_filt$genome)

4.1 Calculate diversities

4.1.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>% 
  remove_rownames() %>% 
  column_to_rownames(var = "genome") %>% 
      filter(rowSums(. != 0, na.rm = TRUE) > 0) %>% 
  rownames_to_column(., "genome")

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt1 %>%
  filter(genome %in% rownames(dist)) %>% 
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample))%>%
  full_join(functional, by = join_by(sample == sample))

treatment_colors <- c('#008080', "#d57d2c")

4.1.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt1 %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

save(beta_q0n, 
     beta_q1n, 
     beta_q1p, 
     beta_q1f, 
     file = "data/beta_27022025.Rdata")

4.2 Does captivity change the microbial community?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") ) 

4.2.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") ) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(7.7073)  ( log )
Formula: richness ~ time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   341.1    347.3   -166.5    333.1       31 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.88922 -0.29614  0.04658  0.45758  1.21929 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.1706   0.413   
Number of obs: 35, groups:  individual, 18

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)               4.1983     0.1330  31.564  < 2e-16 ***
time_point1_Acclimation  -0.4405     0.1344  -3.279  0.00104 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
tm_pnt1_Acc -0.456

Neutral

Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  258.7499 264.7359 -125.3749

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    8.169818 7.178627

Fixed effects:  neutral ~ time_point 
                            Value Std.Error DF   t-value p-value
(Intercept)              36.65610  2.563403 17 14.299780       0
time_point1_Acclimation -16.81743  2.447304 16 -6.871818       0
 Correlation: 
                        (Intr)
time_point1_Acclimation -0.456

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.52249470 -0.58253897 -0.03654165  0.58990386  1.40500351 

Number of Observations: 35
Number of Groups: 18 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  129.7532 135.7392 -60.8766

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.6806603 1.248658

Fixed effects:  phylogenetic ~ time_point 
                           Value Std.Error DF   t-value p-value
(Intercept)             5.221589 0.3351985 17 15.577603  0.0000
time_point1_Acclimation 0.135310 0.4236755 16  0.319371  0.7536
 Correlation: 
                        (Intr)
time_point1_Acclimation -0.61 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.18114503 -0.47484996  0.08008291  0.37909497  1.82560478 

Number of Observations: 35
Number of Groups: 18 

Functional

Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_funct)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -52.38562 -46.39959 30.19281

Random effects:
 Formula: ~1 | individual
        (Intercept)   Residual
StdDev:  0.03084103 0.08367718

Fixed effects:  functional ~ time_point 
                             Value  Std.Error DF  t-value p-value
(Intercept)              1.5154110 0.02101988 17 72.09416  0.0000
time_point1_Acclimation -0.0307069 0.02834791 16 -1.08322  0.2948
 Correlation: 
                        (Intr)
time_point1_Acclimation -0.653

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.20739889 -0.40096259  0.05458715  0.62964551  1.42654682 

Number of Observations: 35
Number of Groups: 18 

4.2.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07277 0.072771 7.2819    999  0.012 *
Residuals 33 0.32979 0.009993                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.5295984	0.05492497	1.917862	0.001
Residual	33	9.1126169	0.94507503	NA	NA
Total	34	9.6422153	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02764 0.027637 2.2078    999  0.128
Residuals 33 0.41309 0.012518                     

adonis2(neutral ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.7470135	0.08504426	3.067318	0.001
Residual	33	8.0368067	0.91495574	NA	NA
Total	34	8.7838201	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.04119 0.041187 2.897    999  0.089 .
Residuals 33 0.46917 0.014217                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2075719	0.1133094	4.217042	0.001
Residual	33	1.6243310	0.8866906	NA	NA
Total	34	1.8319029	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00017 0.0001696 0.0078    999  0.931
Residuals 33 0.71884 0.0217831                     

adonis2(func ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0009227428	0.0009988421	0.03299475	0.615
Residual	33	0.9228897312	0.9990011579	NA	NA
Total	34	0.9238124740	1.0000000000	NA	NA

4.2.3 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

wild_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "0_Wild") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_wild= all(c_across(all_of(wild_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_wild = mean(c_across(all_of(wild_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_wild==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_wild ~ "wild",
      !all_zeros_acclimation & all_zeros_wild ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_wild ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_wild)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 6 × 2
# Rowwise: 
  phylum            sample_count
  <chr>                    <int>
1 p__Pseudomonadota            5
2 p__Bacillota                 3
3 p__Bacillota_A               2
4 p__Bacteroidota              2
5 p__Actinomycetota            1
6 p__Chlamydiota               1

Wild

structural_zeros %>% 
  filter(present=="wild") %>% 
  select(1,5:10)

# A tibble: 11 × 7
# Rowwise: 
   genome                phylum           class           order              family               genus                   species
   <chr>                 <chr>            <chr>           <chr>              <chr>                <chr>                   <chr>  
 1 AH1_2nd_19:bin_000031 p__Bacillota_A   c__Clostridia   o__Lachnospirales  f__Lachnospiraceae   g__Lachnotalea          s__    
 2 AH1_2nd_7:bin_000023  p__Bacillota_A   c__Clostridia   o__Lachnospirales  f__                  g__                     s__    
 3 AH1_2nd_20:bin_000025 p__Bacillota_A   c__Clostridia   o__Lachnospirales  f__CAG-274           g__                     s__    
 4 AH1_2nd_6:bin_000024  p__Spirochaetota c__Brevinematia o__Brevinematales  f__Brevinemataceae   g__Brevinema            s__    
 5 AH1_2nd_10:bin_000097 p__Bacteroidota  c__Bacteroidia  o__Bacteroidales   f__Tannerellaceae    g__Parabacteroides      s__    
 6 LI1_2nd_7:bin_000020  p__Bacillota_B   c__Peptococcia  o__Peptococcales   f__Peptococcaceae    g__                     s__    
 7 LI1_2nd_7:bin_000001  p__Bacillota_A   c__Clostridia   o__Clostridiales   f__Clostridiaceae    g__Sarcina              s__    
 8 LI1_2nd_7:bin_000083  p__Bacillota_A   c__Clostridia   o__Oscillospirales f__Oscillospiraceae  g__Pseudoflavonifractor s__    
 9 AH1_2nd_7:bin_000055  p__Bacillota     c__Bacilli      o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma           s__    
10 LI1_2nd_1:bin_000052  p__Bacillota_A   c__Clostridia   o__TANB77          f__CAG-508           g__RGIG8482             s__    
11 AH1_2nd_6:bin_000060  p__Bacillota_A   c__Clostridia   o__Oscillospirales f__Oscillospiraceae  g__Intestinimonas       s__    

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 35 × 3
   trait    p_value p_adjust
   <chr>      <dbl>    <dbl>
 1 B0214 0.000261   0.00340 
 2 B0219 0.00164    0.0111  
 3 B0302 0.000198   0.00284 
 4 B0703 0.00306    0.0182  
 5 B0709 0.000103   0.00284 
 6 B0712 0.000551   0.00561 
 7 B0801 0.000198   0.00284 
 8 B1012 0.00136    0.0102  
 9 B1028 0.00000837 0.000599
10 D0102 0.00131    0.0102  
# ℹ 25 more rows

4.2.4 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, 
     file = "data/ancombc_wild_accl.Rdata")

load("data/ancombc_wild_accl.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
  filter(p_time_point1_Acclimation < 0.05) %>%
  dplyr::arrange(p_time_point1_Acclimation) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point1_Acclimation)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

        phylum sample_count
1  Bacillota_A           16
2 Bacteroidota            9
3  Bacillota_C            1

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  select(phylum, genus, species)

         phylum                genus                      species
1   Bacillota_A           MGBC140009                             
2   Bacillota_A                                                  
3   Bacillota_A                                                  
4   Bacillota_A     Negativibacillus                             
5  Bacteroidota      Parabacteroides                             
6  Bacteroidota          Bacteroides                             
7   Bacillota_A           Hungatella                             
8  Bacteroidota          Bacteroides                             
9  Bacteroidota      Parabacteroides                             
10 Bacteroidota      Parabacteroides                             
11 Bacteroidota          Bacteroides                             
12  Bacillota_A                                                  
13 Bacteroidota          Bacteroides                             
14  Bacillota_A                                                  
15  Bacillota_A           Copromonas                             
16  Bacillota_A        Enterocloster                             
17  Bacillota_A      Velocimicrobium                             
18  Bacillota_A               CAG-95                             
19  Bacillota_C                                                  
20 Bacteroidota          Bacteroides                             
21  Bacillota_A                                                  
22  Bacillota_A           Copromonas                             
23 Bacteroidota          Bacteroides Bacteroides thetaiotaomicron
24  Bacillota_A Pseudoflavonifractor                             
25  Bacillota_A        Enterocloster                             
26  Bacillota_A             JALFVM01                             

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))

                  genus sample_count
1                                  6
2           Bacteroides            6
3       Parabacteroides            3
4            Copromonas            2
5         Enterocloster            2
6                CAG-95            1
7            Hungatella            1
8              JALFVM01            1
9            MGBC140009            1
10     Negativibacillus            1
11 Pseudoflavonifractor            1
12      Velocimicrobium            1

4.2.5 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Wild"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.2.5.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.331 0.0320
2 Wild        0.345 0.0234

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94347, p-value = 0.07144

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 114, p-value = 0.2066
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Wild)%>% 
  mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.2.5.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.329 0.0320
2 Wild        0.346 0.0195

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94492, p-value = 0.07913

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 100, p-value = 0.08313
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Wild	Function	Difference	treatment
B03	0.1097356	0.1256091	Amino acid derivative biosynthesis	-0.0158735	Wild
D03	0.2921075	0.3442827	Sugar degradation	-0.0521752	Wild

4.3 Does the antibiotic administration change the microbial community?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) 

4.3.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(neutral ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.0665)  ( log )
Formula: neutral ~ time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   223.0    228.8   -107.5    215.0       28 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4702 -0.6482 -0.1777  0.6405  2.6346 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.07573  0.2752  
Number of obs: 32, groups:  individual, 18

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)               2.9069     0.1546  18.800  < 2e-16 ***
time_point2_Antibiotics  -0.9393     0.2131  -4.407 1.05e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
tm_pnt2_Ant -0.538

Neutral

Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  224.1303 229.7351 -108.0652

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.198228 7.479481

Fixed effects:  neutral ~ time_point 
                            Value Std.Error DF   t-value p-value
(Intercept)              19.11841  1.971629 17  9.696759   0e+00
time_point2_Antibiotics -11.46314  2.675428 13 -4.284600   9e-04
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.63 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.86893024 -0.51778639 -0.08383536  0.57135437  1.94021414 

Number of Observations: 32
Number of Groups: 18 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  137.9324 143.5371 -64.96618

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    1.019275  1.66934

Fixed effects:  phylogenetic ~ time_point 
                            Value Std.Error DF   t-value p-value
(Intercept)              5.400123 0.4734177 17 11.406677   0.000
time_point2_Antibiotics -0.791476 0.6015713 13 -1.315681   0.211
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.587

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.65991143 -0.50025273 -0.03029114  0.42477766  1.61395665 

Number of Observations: 32
Number of Groups: 18 

Functional

Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_funct)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -15.05258 -9.447793 11.52629

Random effects:
 Formula: ~1 | individual
         (Intercept)  Residual
StdDev: 1.585329e-06 0.1502428

Fixed effects:  functional ~ time_point 
                             Value  Std.Error DF  t-value p-value
(Intercept)              1.4848858 0.03643923 17 40.74965   0.000
time_point2_Antibiotics -0.0341744 0.05322290 13 -0.64210   0.532
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.685

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.2441309 -0.5590927  0.2912519  0.6141833  2.0349226 

Number of Observations: 32
Number of Groups: 18 

4.3.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.013687 0.0136873 2.0494    999  0.162
Residuals 30 0.200356 0.0066785                     

adonis2(richness ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.580367	0.1313608	4.536779	0.001
Residual	30	10.450370	0.8686392	NA	NA
Total	31	12.030738	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03491 0.034906 2.5203    999  0.136
Residuals 30 0.41549 0.013850                     

adonis2(neutral ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.78013	0.1607184	5.744856	0.001
Residual	30	9.29595	0.8392816	NA	NA
Total	31	11.07608	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.31336 0.313365 16.503    999  0.002 **
Residuals 30 0.56964 0.018988                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.572915	0.2890287	12.19579	0.001
Residual	30	3.869157	0.7109713	NA	NA
Total	31	5.442071	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18493 0.184935 4.9256    999  0.037 *
Residuals 30 1.12638 0.037546                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(func ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.120331	0.3613445	16.97368	0.001
Residual	30	1.980120	0.6386555	NA	NA
Total	31	3.100451	1.0000000	NA	NA

4.3.3 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "2_Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               7
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

4.3.4 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, 
     file = "data/ancombc_antib_accl.Rdata")

load("data/ancombc_antib_accl.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
  filter(p_time_point2_Antibiotics < 0.05) %>%
  dplyr::arrange(p_time_point2_Antibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point2_Antibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

            phylum sample_count
1     Bacteroidota           14
2      Bacillota_A            4
3        Bacillota            2
4 Campylobacterota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  select(phylum, genus, species)

             phylum           genus species
1  Campylobacterota  Helicobacter_J        
2       Bacillota_A                        
3      Bacteroidota     Bacteroides        
4      Bacteroidota     Odoribacter        
5      Bacteroidota     Bacteroides        
6         Bacillota   Coprobacillus        
7         Bacillota                        
8      Bacteroidota     Phocaeicola        
9      Bacteroidota     Bacteroides        
10     Bacteroidota     Phocaeicola        
11     Bacteroidota     Odoribacter        
12     Bacteroidota     Bacteroides        
13     Bacteroidota     Bacteroides        
14     Bacteroidota Parabacteroides        
15     Bacteroidota                        
16     Bacteroidota Parabacteroides        
17      Bacillota_A                        
18     Bacteroidota       Alistipes        
19      Bacillota_A                        
20      Bacillota_A                        
21     Bacteroidota     Odoribacter        

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))

            genus sample_count
1                            6
2     Bacteroides            5
3     Odoribacter            3
4 Parabacteroides            2
5     Phocaeicola            2
6       Alistipes            1
7   Coprobacillus            1
8  Helicobacter_J            1

4.3.5 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.3.5.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.331 0.0320
2 Antibiotics 0.244 0.134 

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94332, p-value = 0.09304

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.3.5.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.329 0.0320
2 Antibiotics 0.255 0.126 

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95737, p-value = 0.2324

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	treatment
B06	0.68110280	0.47325490	Organic anion biosynthesis	0.20784790	Acclimation
B02	0.59963040	0.41562100	Amino acid biosynthesis	0.18400940	Acclimation
D02	0.38587770	0.20636760	Polysaccharide degradation	0.17951010	Acclimation
S03	0.27103471	0.09357224	Spore	0.17746247	Acclimation
B01	0.84215140	0.69078910	Nucleic acid biosynthesis	0.15136230	Acclimation
B07	0.44558700	0.30432910	Vitamin biosynthesis	0.14125790	Acclimation
B08	0.44596150	0.32044710	Aromatic compound biosynthesis	0.12551440	Acclimation
D09	0.25533360	0.13438350	Antibiotic degradation	0.12095010	Acclimation
B04	0.54457570	0.42921780	SCFA biosynthesis	0.11535790	Acclimation
D03	0.29210750	0.20698550	Sugar degradation	0.08512200	Acclimation
D05	0.28856080	0.22258820	Amino acid degradation	0.06597260	Acclimation
B10	0.05947806	0.03621687	Antibiotic biosynthesis	0.02326119	Acclimation

4.4 Does antibiotic administration remove the differences found in the two populations?

4.4.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point == "2_Antibiotics" )

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(time_point == "2_Antibiotics" ) %>% 
  ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
      facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.4.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point == "2_Antibiotics") %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point == "2_Antibiotics")

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.015377 0.0153774 6.8934    999  0.011 *
Residuals 21 0.046845 0.0022307                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.361743	0.1533845	3.80465	0.001
Residual	21	7.516224	0.8466155	NA	NA
Total	22	8.877966	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030649 0.0306488 3.8694    999  0.069 .
Residuals 21 0.166339 0.0079209                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.787698	0.208772	5.541022	0.001
Residual	21	6.775221	0.791228	NA	NA
Total	22	8.562919	1.000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012165 0.012165 0.9979    999  0.354
Residuals 21 0.256012 0.012191                     

adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.8970751	0.1890846	4.896659	0.001
Residual	21	3.8472307	0.8109154	NA	NA
Total	22	4.7443057	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01804 0.018036 0.4397    999  0.517
Residuals 21 0.86147 0.041023                     

adonis2(func ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0184103	0.009818749	0.2082384	0.678
Residual	21	1.8566048	0.990181251	NA	NA
Total	22	1.8750151	1.000000000	NA	NA

4.5 Are the microbial communities similar in both donor samples?

samples_to_keep <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))%>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2")) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00014 0.0001396 0.0173    999  0.912
Residuals 13 0.10481 0.0080621                     

adonis2(richness ~ time_point,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.06889268	0.03010712	0.4035421	0.981
Residual	13	2.21935895	0.96989288	NA	NA
Total	14	2.28825162	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000655 0.0006546 0.0514    999  0.816
Residuals 13 0.165409 0.0127237                     

adonis2(neutral ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.09294671	0.03882177	0.5250671	0.818
Residual	13	2.30124351	0.96117823	NA	NA
Total	14	2.39419022	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003691 0.0036912 0.3071    999  0.673
Residuals 13 0.156266 0.0120205                     

adonis2(phylo ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.009773631	0.01921019	0.2546239	0.841
Residual	13	0.498999565	0.98078981	NA	NA
Total	14	0.508773196	1.00000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.013592 0.0135920 1.5644    999  0.195
Residuals 13 0.112945 0.0086881                     

adonis2(func ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	-0.003158048	-0.0167584	-0.2142684	0.821
Residual	13	0.191603701	1.0167584	NA	NA
Total	14	0.188445652	1.0000000	NA	NA

4.6 Does the donor sample maintain the microbial community found in acclimation?

4.6.1 Alpha diversity

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
    TRUE ~ treatment
  ))
samples_to_keep <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.6.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00836 0.0083600 1.4409    999  0.255
Residuals 22 0.12764 0.0058019                     

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2657351	0.06396199	1.503319	0.096
Residual	22	3.8888430	0.93603801	NA	NA
Total	23	4.1545781	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000661 0.0006614 0.0736    999  0.793
Residuals 22 0.197804 0.0089911                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3164816	0.07596593	1.808646	0.07
Residual	22	3.8496178	0.92403407	NA	NA
Total	23	4.1660995	1.00000000	NA	NA

neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
 # scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00204 0.0020405 0.2622    999   0.71
Residuals 22 0.17118 0.0077807                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0412056	0.056244	1.31111	0.268
Residual	22	0.6914166	0.943756	NA	NA
Total	23	0.7326222	1.000000	NA	NA

phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000463 0.0004627 0.0465    999  0.839
Residuals 22 0.219125 0.0099602                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.08931756	0.2172817	6.107175	0.056
Residual	22	0.32175047	0.7827183	NA	NA
Total	23	0.41106803	1.0000000	NA	NA

func %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

4.7 Is the donor sample microbiota different to recipient microbiota?

4.7.1 Alpha diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.7.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.12795 0.127946 13.179    999  0.001 ***
Residuals 20 0.19416 0.009708                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.777815	0.2760196	7.625059	0.001
Residual	20	4.663086	0.7239804	NA	NA
Total	21	6.440902	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.071502 0.071502 5.4967    999  0.031 *
Residuals 20 0.260166 0.013008                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	2.112572	0.3210813	9.458609	0.001
Residual	20	4.466983	0.6789187	NA	NA
Total	21	6.579556	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04731 0.047307 2.6157    999  0.115
Residuals 20 0.36172 0.018086                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3778248	0.239656	6.303884	0.001
Residual	20	1.1987049	0.760344	NA	NA
Total	21	1.5765298	1.000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000167 0.0001673 0.0139    999  0.901
Residuals 20 0.241060 0.0120530                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.07837549	0.1797296	4.382204	0.088
Residual	20	0.35769892	0.8202704	NA	NA
Total	21	0.43607441	1.0000000	NA	NA

4.8 Does FMT change the microbial community over time?

4.8.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1","Post-FMT2") ) 

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type=="Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2" )) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Model_rich <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Model_rich)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.3876)  ( log )
Formula: richness ~ treatment + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   171.0    174.3    -81.5    163.0       13 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.73495 -0.71265 -0.07086  0.40744  1.88240 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups:  individual, 9

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.9343     0.1759  22.369   <2e-16 ***
treatmentPost-FMT2   0.4052     0.2402   1.687   0.0917 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trtmnP-FMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  128.7946 131.6268 -60.3973

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.497472 11.25384

Fixed effects:  neutral ~ treatment 
                      Value Std.Error DF  t-value p-value
(Intercept)        24.65947  4.164756  8 5.920988  0.0004
treatmentPost-FMT2 14.87494  5.482531  7 2.713151  0.0301
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.7  

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.76462323 -0.55701822 -0.04763166  0.55267588  1.30333436 

Number of Observations: 17
Number of Groups: 9 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  56.17121 59.00341 -24.0856

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.6979711 0.8383981

Fixed effects:  phylogenetic ~ treatment 
                      Value Std.Error DF   t-value p-value
(Intercept)        4.163606 0.3820859  8 10.897043  0.0000
treatmentPost-FMT2 0.916226 0.4122640  7  2.222426  0.0617
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.583

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.1362905 -0.5779277 -0.1466831  0.4812394  2.3357705 

Number of Observations: 17
Number of Groups: 9 

Functional

Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_funct)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -32.72222 -29.89002 20.36111

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:  0.06160077 0.0305115

Fixed effects:  functional ~ treatment 
                       Value  Std.Error DF  t-value p-value
(Intercept)        1.4877804 0.02341751  8 63.53281  0.0000
treatmentPost-FMT2 0.0320561 0.01517203  7  2.11284  0.0725
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.357

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.07361856 -0.66739366 -0.04576292  0.53118788  1.61981397 

Number of Observations: 17
Number of Groups: 9 

4.8.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.017610 0.017610 2.8225    999  0.102
Residuals 15 0.093585 0.006239                     

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3560084	0.07968734	1.298809	0.00390625
Residual	15	4.1115570	0.92031266	NA	NA
Total	16	4.4675655	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.009828 0.0098278 0.7113    999  0.436
Residuals 15 0.207263 0.0138175                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3149954	0.08110541	1.323962	0.00390625
Residual	15	3.5687823	0.91889459	NA	NA
Total	16	3.8837776	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.010955 0.010955 0.6602    999  0.545
Residuals 15 0.248892 0.016593                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.06391536	0.09449338	1.565312	0.0703125
Residual	15	0.61248501	0.90550662	NA	NA
Total	16	0.67640037	1.00000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00151 0.0015101 0.1075    999  0.727
Residuals 15 0.21063 0.0140417                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.02861781	0.07980286	1.300855	0.3203125
Residual	15	0.32998847	0.92019714	NA	NA
Total	16	0.35860627	1.00000000	NA	NA

4.9 Do FMT change the microbial community compared to antibiotics and acclimation?

4.9.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","Post-FMT1") ) 

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "Post-FMT1")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
#  scale_color_manual(values=treatment_colors)+
#  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness: Problems to calculate

Model_rich <- glmer.nb(richness ~ treatment + (1|individual), data = alpha_div_meta)
#summary(Model_rich)
emmeans(Model_rich, pairwise ~ treatment)

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_neutral)
emmeans(Model_neutral, pairwise ~ treatment)

$emmeans
 treatment   emmean   SE df lower.CL upper.CL
 Acclimation  17.31 3.42  8    9.431     25.2
 Antibiotics   8.44 3.86  8   -0.451     17.3
 Post-FMT1    24.60 3.62  8   16.256     32.9

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate   SE df t.ratio p.value
 Acclimation - Antibiotics     8.87 4.74 13   1.869  0.1868
 Acclimation - (Post-FMT1)    -7.29 4.55 13  -1.601  0.2800
 Antibiotics - (Post-FMT1)   -16.16 4.85 13  -3.333  0.0139

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_phylo)
emmeans(Model_phylo, pairwise ~ treatment)

$emmeans
 treatment   emmean    SE df lower.CL upper.CL
 Acclimation   5.50 0.505  8     4.34     6.67
 Antibiotics   4.30 0.571  8     2.98     5.62
 Post-FMT1     4.15 0.535  8     2.91     5.38

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate    SE df t.ratio p.value
 Acclimation - Antibiotics    1.203 0.716 13   1.680  0.2496
 Acclimation - (Post-FMT1)    1.357 0.688 13   1.973  0.1583
 Antibiotics - (Post-FMT1)    0.154 0.733 13   0.210  0.9760

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

Functional

Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_funct)
emmeans(Model_funct, pairwise ~ treatment)

$emmeans
 treatment   emmean     SE df lower.CL upper.CL
 Acclimation   1.49 0.0437  8     1.39     1.60
 Antibiotics   1.47 0.0495  8     1.35     1.58
 Post-FMT1     1.49 0.0463  8     1.38     1.60

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate     SE df t.ratio p.value
 Acclimation - Antibiotics  0.02530 0.0660 13   0.383  0.9227
 Acclimation - (Post-FMT1)  0.00618 0.0637 13   0.097  0.9948
 Antibiotics - (Post-FMT1) -0.01913 0.0678 13  -0.282  0.9572

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

4.9.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq    F N.Perm Pr(>F)
Groups     2 0.01582 0.0079101 1.05    999  0.383
Residuals 21 0.15821 0.0075336                   

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.840841	0.2018508	2.655435	0.001
Residual	21	7.278967	0.7981492	NA	NA
Total	23	9.119808	1.0000000	NA	NA

pairwise.adonis(richness, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8287871 2.293722 0.1407734   0.002      0.006   *
2   Acclimation vs Post-FMT1  1 0.9572484 2.939042 0.1638350   0.001      0.003   *
3   Antibiotics vs Post-FMT1  1 0.9764237 2.751191 0.1746656   0.001      0.003   *

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01796 0.0089778 0.5959    999  0.592
Residuals 21 0.31640 0.0150667                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	2.357142	0.2699241	3.882065	0.001
Residual	21	6.375470	0.7300759	NA	NA
Total	23	8.732613	1.0000000	NA	NA

pairwise.adonis(neutral, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8814932 2.747044 0.1640316   0.001      0.003   *
2   Acclimation vs Post-FMT1  1 1.3133104 4.634354 0.2360329   0.001      0.003   *
3   Antibiotics vs Post-FMT1  1 1.3427497 4.355528 0.2509591   0.002      0.006   *

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.15410 0.077052 2.6434    999   0.08 .
Residuals 21 0.61214 0.029149                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.327076	0.3647564	6.029092	0.001
Residual	21	2.311178	0.6352436	NA	NA
Total	23	3.638254	1.0000000	NA	NA

pairwise.adonis(phylo, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.6885926 5.124832 0.2679674   0.005      0.015   .
2   Acclimation vs Post-FMT1  1 0.2548027 3.292748 0.1800029   0.045      0.135    
3   Antibiotics vs Post-FMT1  1 1.1000473 9.048069 0.4103792   0.002      0.006   *

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.23408 0.117038 4.3947    999  0.024 *
Residuals 21 0.55927 0.026632                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.636530	0.5730275	14.09175	0.001
Residual	21	1.219406	0.4269725	NA	NA
Total	23	2.855936	1.0000000	NA	NA

pairwise.adonis(func, subset_meta$treatment, perm = 999)

                       pairs Df  SumsOfSqs   F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 1.11236925 14.811924 0.5140901   0.003      0.009   *
2   Acclimation vs Post-FMT1  1 0.05914673  2.630593 0.1492061   0.166      0.498    
3   Antibiotics vs Post-FMT1  1 1.36488774 16.896102 0.5651607   0.001      0.003   *

4.10 Is the gut microbiota similar to the donor after FMT?

4.10.1 Beta diversity

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  select(sample) %>% 
  pull()

subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.11100 0.111002 11.752    999  0.001 ***
Residuals 28 0.26447 0.009445                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.184578	0.1548451	5.13002	0.001
Residual	28	6.465508	0.8451549	NA	NA
Total	29	7.650086	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04408 0.044082 3.1744    999  0.092 .
Residuals 28 0.38883 0.013887                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.147077	0.1608783	5.368223	0.001
Residual	28	5.983012	0.8391217	NA	NA
Total	29	7.130089	1.0000000	NA	NA

neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
 # scale_color_manual(values = treatment_colors) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00007 0.0000689 0.0044    999  0.948
Residuals 28 0.43637 0.0155847                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.1298187	0.1018776	3.17615	0.001
Residual	28	1.1444432	0.8981224	NA	NA
Total	29	1.2742619	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00121 0.0012109 0.1039    999   0.78
Residuals 28 0.32638 0.0116565                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.01515762	0.02696066	0.7758149	0.481
Residual	28	0.54705472	0.97303934	NA	NA
Total	29	0.56221234	1.00000000	NA	NA

4.10.2 Structural zeros

FMT_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment == "FMT") %>% 
  dplyr::select(sample) %>%
  pull()

Transplant_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment =="Transplant") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>% 
   mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>% 
   mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>% 
   mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
      !all_zeros_FMT & all_zeros_Transplant ~ "FMT",
      !all_zeros_FMT & !all_zeros_Transplant ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Structural zeros in each group

fmt <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

fmt_f <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(family, name = "FMT") %>%
  arrange(desc(FMT)) 

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>%
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))

  Transplant FMT
1         52  55

Phyla to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()

phylum	Transplant	FMT
p__Bacillota_A	20	21
p__Bacillota	16	10
p__Pseudomonadota	6	11
p__Bacteroidota	3	6
p__Desulfobacterota	2	2
p__Bacillota_B	1	0
p__Chlamydiota	1	0
p__Cyanobacteriota	1	0
p__Fusobacteriota	1	0
p__Verrucomicrobiota	1	2
p__Actinomycetota	0	1
p__Bacillota_C	0	1
p__Campylobacterota	0	1

Families to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(family, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt_f, by="family" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()

family	Transplant	FMT
f__Lachnospiraceae	7	9
f__Erysipelotrichaceae	6	5
f__UBA660	6	0
f__Enterobacteriaceae	5	2
f__CAG-508	3	0
f__Ruminococcaceae	3	3
f__Anaerovoracaceae	2	0
f__Coprobacillaceae	2	0
f__Desulfovibrionaceae	2	2
f__Oscillospiraceae	2	1
f__Tannerellaceae	2	1
f__UBA1242	2	0
f__	1	3
f__Akkermansiaceae	1	0
f__CAG-239	1	2
f__Chlamydiaceae	1	0
f__Enterococcaceae	1	3
f__Fusobacteriaceae	1	0
f__Gastranaerophilaceae	1	0
f__Marinifilaceae	1	0
f__Mycoplasmoidaceae	1	1
f__Peptococcaceae	1	0
f__Anaerotignaceae	0	2
f__Bacteroidaceae	0	2
f__LL51	0	2
f__UBA3700	0	2
f__Acutalibacteraceae	0	1
f__Arcobacteraceae	0	1
f__CAG-274	0	1
f__CALVMC01	0	1
f__Devosiaceae	0	1
f__Mycobacteriaceae	0	1
f__Pumilibacteraceae	0	1
f__RUG11792	0	1
f__Rhizobiaceae	0	1
f__Rikenellaceae	0	1
f__Sphingobacteriaceae	0	1
f__Streptococcaceae	0	1
f__UBA1997	0	1
f__UBA3830	0	1
f__Weeksellaceae	0	1

4.10.2.1 Functional capacities of the structural zeros

#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  filter(genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
  column_to_rownames(., "genome") 
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Elevation")

4.10.3 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  column_to_rownames("sample") %>%
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "treatment",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, file = "data/ancombc_FMT_transplant.Rdata")

load("data/ancombc_FMT_transplant.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
  filter(p_treatmentTransplant < 0.05) %>%
  dplyr::arrange(p_treatmentTransplant) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_treatmentTransplant)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

Differential abundance MAGs in each group

fmt_count <- ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant<0)  %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0)))

             phylum Transplant FMT
1      Bacteroidota         13   5
2       Bacillota_A          4  17
3         Bacillota          3   1
4  Campylobacterota          1   1
5  Desulfobacterota          1   2
6     Spirochaetota          1   0
7    Pseudomonadota          0   5
8   Cyanobacteriota          0   2
9       Bacillota_B          0   1
10      Bacillota_C          0   1

ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0))) %>% 
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))

  Transplant FMT
1         23  35

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

4.10.4 Differences in the functional capacities

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.10.4.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.354 0.0255
2 Transplant 0.351 0.0426

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94871, p-value = 0.1561

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 

means_gift <- element_gift_filt %>% 
  select(-treatment) %>% 
  pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>% 
  left_join(sample_metadata, by=join_by(sample==sample)) %>% 
  group_by(treatment, Elements) %>%
  summarise(mean=mean(abundance))

log_fold <- means_gift %>%
  group_by(Elements) %>%
  summarise(
    logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
    ) %>% 
  left_join(., difference_table, by="Elements")

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

log_fold %>%
  filter(Elements!="D0611") %>% 
  ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Log_fold")+
  labs(fill="Treatment")

4.10.4.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.349 0.0221
2 Transplant 0.354 0.0375

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.87207, p-value = 0.001864

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 120, p-value = 0.7106
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))